New algorithm for classical gauge theory simulations in an expanding box 
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We propose a new algorithm for classical statistical simulations in scalar and gauge theories 
undergoing a one dimensional expansion, which allows simulations to study boxes of larger transverse 
extent and to continue for longer times, without losing lattice resolution in the expanding direction. 
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I. INTRODUCTION 

The initial condition of heavy ion collisions of asymp- 
totically large nuclei at asymptotically large energies is 
described by the color glass condensate framework [3, 0] , 
in which the relevant degrees of freedom right after the 
collision are nearly boost invariant large-amplitude gauge 
fields, or "glasma" fields [|[. The subsequent evolution 
towards local thermal equilibrium is of phcnomenological 
interest and has attracted a lot of theoretical attention. 

The large amplitudes of the fields admit a classical sta- 
tistical treatment, and the evolution of the fields at early 
times can be followed in classical statistical lattice simu- 
lations. In the case of exact boost invariance, it is enough 
to perform simulations in 2 + ID However, fluctu- 

ations are not boost invariant, and some fluctuations are 
unstable to exponential growth Q. Therefore the study 
of this system, and particularly of the growth and fate of 
these fluctuations, requires classical lattice gauge theory 
studies in longitudinally expanding 3+1 dimensions. 

In the simulations, one discretizes the gauge fields, and 
solves the classical Yang-Mills equations for time evolu- 
tion. The approximate statistical boost invariance of the 
system makes it most convenient to discretize the co- 
moving coordinates r = \/t 2 — z 2 and rj = atanh(z/f). 
In these coordinates the physical lattice spacing in the 
longitudinal direction a n = tAtj grows linearly as a func- 
tion of proper time, while the transverse lattice spacing 
aj_ stays constant. This, however, introduces the bottle- 
neck numerical challenge of these simulations: the lon- 
gitudinal spacing tAtj cannot be much larger than the 
transverse spacing even at the end of the simulation Tf ; 
but this requires that one start the simulation with an 
extremely fine lattice in the ^-direction, and hence with 
a fantastically large aspect ratio N v /N± to have enough 
dynamical range in proper time. For a given number of 
lattice points available, this inevitably restricts simula- 
tions to boxes with small physical transverse size, making 
it very difficult to accommodate all the physical scales 
inside the lattice. Indeed, the state of the art simula- 
tions in expanding lattices are being performed on lat- 
tices with e.g., Nj_ x N v = 32 2 x 1024 [|, which is to 
be compared with simulations in static boxes reaching 
TV 3 = 256 3 This is especially problematic because 
the classical field evolution becomes a multi-scale prob- 
lem; there is the scale Q s associated with the structure 
of the initial conditions, and there is a screening scale m 



which at late times is parametrically m ~ Q s (Q s t) -1 / 2 . 
The transverse lattice spacing needs to be fine enough 
to resolve the scale Q s , Q s a± < 1, but the transverse 
size should be enough to contain the scale to, mL± > 1 
(with L± = N±a±_); otherwise important features of the 
dynamics may be missed. 

The purpose of this short note is to propose a new algo- 
rithm to overcome this problem, and to facilitate simula- 
tions with arbitrary dynamical range in time while having 
large N± . The algorithm works by cropping the lattice in 
half in the ?y-direction whenever the ratio a^/a^ = £ has 
become too large. This procedure reduces the number 
of degrees of freedom to half, which are subsequently re- 
covered by a mesh refinement in the 77-direction, halving 
the physical lattice spacing a v . In non-abelian gauge the- 
ories, the mesh refinement procedure generates Gauss's 
law violations of the order of a 2 F, which arc then elim- 
inated in the third step of the algorithm by projecting 
the new field configuration to the physical manifold. 

In Section [Til we describe our algorithm in a simpler 
scalar theory. For the scalar theory, the implementation 
is extremely simple and should be easily incorporated 
into any lattice practitioner's code without effort. We 
then move to non-abelian gauge theory and consider the 
subtleties that arise in that case. 



II. SCALAR FIELD THEORY 

As a warmup, let us consider a lattice scalar field the- 
ory discretized on a co-moving lattice with a lattice spac- 
ing a± in the transverse direction and A77 in the rapidity 
direction. The extent of the lattice is Nj_ x (N v + 1); the 
field <f>fj,a and its conjugate momentum 7r^n live on the 
lattice sites x = (77, hi, n.2), labeled by the indices 77 = 
{0,...,^} andn= ({0, . . . , N ± - 1}, {0, . . . , N± - 1}) 
in rapidity and transverse directions, respectively. Their 
evolution in proper time is given by equations of motion 
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for a generic potential V. The operators and d 2 are 
some implementation of lattice Laplacians, in the sim- 
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plcst case just symmetric differences 
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where the summation goes over the transverse directions. 
One way to obtain this form is by deriving it from a 
Hamiltonian written in terms of the lattice fields, 
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where a 2 L rA?7 ^2 — J d x±Tdr) and the three terms in 
parenthesis are 7r 2 , t~ 2 (<9 ?) 0) 2 , and (Vj_</>) 2 respectively. 

In practical simulations, one typically imposes periodic 
boundary conditions (BC) to minimize the finite volume 
effects and for the transverse directions this is our choice. 
But in the rapidity direction we find it beneficial to im- 
pose Neumann BC, that is, we discard the term in Eq. ([5]) 
involving (0jy,,+i=o,n - </>iv„,n) 2 or equivalently set 
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Like periodic BC, Neumann BC conserve the energy 1 . 
And we will find that the cutting of the lattice is better 
implemented with Neumann BC. 

Our algorithm works as follows: Whenever the ratio of 
lattice spacings £ = a v /aj_ reaches a fiducial value £ c < 1 
we do the following. First, we divide the lattice into three 
rapidity regions {0, . . . , N bcl - 1}, {N bcl , N bc2 }, and 
{N bc2 + 1 . . . JV ? }, with N bcl = JV„/4 and iV 6c2 = 3N„/4 
as shown in Figure [1] We discard the end caps, and 
impose Neumann BC for the middle region 
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This reduces the length of our lattice by a factor of 2. 

In the second step (carried out at the same time r), 
we refine the mesh in the rapidity direction. Consider a 
new lattice with half the lattice spacing in the rapidity 
Qfj/2, such that 
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= 2(f,-N bcl ) = {Q,...,N v }. 
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FIG. 1: In the first step of the algorithm, the extent of the 
lattice in the rapidity direction is reduced by a factor 2 by 
cropping the ends of the lattice. 



The original lattice field specifies the value of the new 
fields and momenta only at lattice sites with even fj new . 
Thus we need to provide a prescription how to interpo- 
late the fields at the odd lattice sites, and the prescription 
should be such that it keeps the fields as smooth as possi- 
ble, that is, it does not transfer energy from the infrared 
modes to the corners of the Brillouin zone. In the case 
of scalar theory, such a description is trivial and in the 
simplest case one can just linearly interpolate fields in 
the rapidity direction 2 
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In terms of the new 77-variable, the Neumann BC reads 
as in Eq. This completes the algorithm. 

Without using our algorithm, there are lattice spac- 
ing errors which scale as 0{a 2 i _k 2 A _) and as C(a 2 fc 2 ) = 
0(A?7 2 fc 2 T 2 ). The latter error grows with time. By re- 
fining the lattice, we prevent this growth with time. How- 
ever, there are three places where new systematic errors 
arise, though each proves to be manageable. Firstly, 
using Neumann BC physically corresponds to replac- 
ing the end of the lattice with a mirror, leading to un- 
physical interference effects within one correlation length 
(A(t]t) ~ of the boundary. This can be easily ame- 

liorated by performing all measurements only in a fiducial 
volume away from the boundary. Very conservatively, 
one may define expectation values by 
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and use lattices with aspect ratio \ = N^/N^ > 2. 

Secondly, abruptly imposing Neumann BC creates a 
configuration in the new lattice which has a cusp in the 
fields where the new boundary lies. The cusp contains 
unphysical ultraviolet modes which subsequently propa- 
gate to the region where measurements are performed 



1 One could equally well use Dirichlet BC and fix the field value 
at the boundary. However we expect (4?) to shrink with time, 
so this may lead to larger finite size effects. 



2 Done properly, the lattice cutting and interpolation can be per- 
formed in the same computer memory as the fields <j>, it are stored 
for the evolution. 
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and thus cause contamination. However, one can in- 
troduce the new BC adiabatically as follows. Rather 
than abruptly introducing "mirrors" at Nf, c i and Nf, C 2, 
one gradually introduces "partly silvered" mirrors, with 
silvering fraction Ag. Specifically, one multiplies the 

(^N bca +i,A - 0Af bc2 ,n) 2 and (4>N bcl ,A - <?W bcl -i,n) 2 terms 
in Eq. ([5]) by a coefficient [f — Ag(r)], which smoothly 
changes from Ag(r) = a few correlation times before 
the cut is to be made, to Ag(r) = 1 at the time when 
the cut is to be made. When Ag(r) = (zero silver- 
ing) the lattice obeys the usual update rules. When 
Ag(r) = 1 (full silvering), there are Neumann BC's at 
Nbd and Nb C 2- For intermediate values, waves propagat- 
ing towards the Nb c i plane from either side will partially 
reflect as from a half-silvered mirror, with Ag(r) the ex- 
tent of "silvering" of the mirror. Once Ag(r) = 1, equa- 
tions of motion for points between f] = Nb c i and fj = Nb C 2 
inclusive no longer make any reference to points outside 
this range, so the cut is then harmless. 

Thirdly, in performing the interpolation, one system- 
atically underestimates the energy density by an 0(alkZ) 
amount. This is, however, a subleading source of error in 
the simulation compared to the discretization error from 
the transverse directions O(o 2 1 k 2 _) for two reasons. First, 
as long as one keeps the fiducial ratio £ c < 1, the lattice 
spacing in the rj direction is smaller that the transverse 
one. And secondly, due to the longitudinal expansion, 
k v gets redshifted and one generally expects k n < k±. 
Hence, during the whole simulation, up to arbitrarily late 
times, the systematical errors are bounded by those aris- 
ing from the discretization of the transverse direction. 



III. GAUGE THEORY 

Now we apply these ideas to nonabelian gauge theory 
defined on the same anisotropic expanding lattice. The 
Hamiltonian (in A T = gauge) is 3 

h = «i«.E{ E Tr K" 2 ^ 2 J 

77, n 1—77,1,2 
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where the plaquette operator □ is the product of the link 
matrices around an elementary plaquette 
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3 Our normalization is related to that of Q by _Ejj Gre = c a /a T E a 
in the notation of the reference. The usual continuum normal- 



where the link matrices are elements of the group, 
and their canonical momenta E^^ belong to the Lie alge- 
bra of the group. The corresponding equations of motion 
read 
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where ad[C/] a = 2ReTr t a U denotes the adjoint represen- 
tation, where t a are the generators of the lie algebra in 
canonical normalization Trt a t = \S ab , and 



(15) 



The Neumann BC is implemented by setting the links 
and electric fluxes penetrating the boundary to zero 



U\ A = 0, 



E\ A = 0, 
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= 0, El ,=0,(16) 



or equivalently by discarding the terms including these 
links from the Hamiltonian. 

In analogy to the scalar field theory, the first step is 
to cut the lattice in half adiabatically by modifying the 
Hamiltonian. This is done by multiplying the terms in 
Hamiltonian containing plaqucttcs connecting the differ- 
ent rapidity regions (Q^ cl _ lifi and Q^ o2>ft ) by the sil- 
vering function [1 — Ag(r)] which again smoothly goes 
from Ag(r) = to Ag(r) = 1. This modification makes 
sure that the resulting configuration, when the mirror 
becomes fully reflecting, is smooth. 

An extra complication that arises in the gauge theory 
comes from the need to satisfy Gauss' law 
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Ei-Ut ia Ei-e a US_ ia )=0. (17) 



Before cropping the lattice, Gauss' law at the site 
(iVb C 2, n) receives a contribution from 2 fi . Abruptly 
discarding everything to the right of Nb C 2 will discard 
this flux, and Gauss' law will not be satisfied on the 
Nbc2 surface (or the N^d surface). Physically this means 
that there will be "charges" trapped on the surface, rep- 
resenting the flux which previously propagated into the 
mirror. The flux can, however, be forced to go to zero by 
also adiabatically removing the terms in the Hamiltonian 
containing color-electric fields penetrating the mirrors, 
i.e., by multiplying the terms containing {E^ i _ 1 A ) 2 
and (Ex A ) 2 by the function [1 - Ag(r)]. The con- 
tribution of these electric fields to Eq. (fP7|) get multi- 
plied by [1 — Ag(r)], while Eq. (fH| changes from involv- 
ing £fE n N Jar, to involving £[1 - Ag^E^ Ja,,, 
which makes the E field grow correspondingly, preserv- 
ing Gauss' law. But wave evolution will mean that the 
E field naturally evolves to remain about the same size, 
so its contribution to Gauss' law will shrink to zero as 
Ag(r) approaches 1. 
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FIG. 2: The interpolation of the link matrices. The thin dot- 
ted arrows correspond to the links of the original lattice, while 
the thick arrows denote link matrices on the new finer lattice. 
The the color coding (if available) indicates the assignments 
of the old matrices to the new ones. Links in the ^-direction 
starting from odd »7-sites on the the new lattice are set to unit 
matrices (light blue arrows) while the transverse links start- 
ing from odd 77-sites are unitarized averages of the parallel 
transports along the paths shown by the dashed lines. 
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FIG. 3: Interpolation of the transverse electric fields emanat- 
ing from the new lattice sites. The graphical representation 
emphasizes that while under gauge transformations the in- 
fields transform as objects "living" on sites, they are conju- 
gate variables of the gauge links and hence live between sites. 
The transverse i5-field on the new lattice site is the aver- 
age of electric fields on neighboring sites, parallel transported 
along the paths shown by the arrows. The fundamental (anti- 
fundamental) indices are parallel transported along the short 
(long) paths. 



In the second step of the algorithm, the mesh in 
the rapidity direction is refined and fields are interpo- 
lated. We begin by describing the procedure to inter- 
polate the link matrices followed by the electric fields. 
The process of interpolating the gauge field links is il- 
lustrated in Figure HJ The transverse links between pre- 
viously existing lattice sites are left unchanged; that is, 
we choose A = U^ +N j fl for i = 1, 2, where we de- 
note the fields on the refined lattice with tildes, and give 
their coordinates in terms of the refined rapidity vari- 
able ?7 ncw = 2(f) — Nbd). Next consider the two sites 
(n + Nbci,n) and (n + N^i + 1, ft) on the old lattice, 
which become (2n,n) and (2n + 2,n) on the refined lat- 
tice. The new 77-links should obey 

U 2n,{i U 2n+l,A = U n+N bcl ,A ( 18 ) 

so that comparisons between the sites on the refined lat- 
tice agree with those on the unrefined lattice. There is a 
new gauge freedom on the introduced point (2n + l,fi), 
and we will use up that freedom to choose U2 n+1 ft = 1 

the identity. Then = U2 +Nbcl , & . 

Next we interpolate the transverse links of form 
L^2ra+i A , i — 1,2. We take the connection between the 
site (2n + l,n) and the site (2n + l,n + e,) to be the 
re-unitarized average of the connections along the two 
closest paths which use the 77-links and the pre-existing 
transverse links, as shown in Figure [5] 

^2n+l,fi = Pr °jsU(Ar)2 (^2n+l,n^2n+2,fi^2ri t +l,n+e i 

+U?tA A Ul A+ ^) , (19) 
where Projgu^ means projection onto the nearest group 



element, which is just a rescaling in SU(2) and can be 
performed for SU(7V > 2) using the algorithm of Ref. [l(| ■ 

With the links defined, we turn to the electric fields. 
Though Eq. (Ti"4"|) defines the electric field E? fi as an ad- 
joint object transforming at the site (77, ft), the i? Q -fields 
are conjugate variables of the link variables U a and hence 
it is most natural to think of them as "living" on the links; 
the natural objects are in fact E a U a which belong to the 
tangent space of U a rather than the lie algebra elements 
E a . 

We choose the E^-ficlds on the existing transverse links 
to be unchanged, E l 2n h = E^ n+N h . The E^-fields on 
the refined lattice, viewed as living at the centerpoints 
of their links, lie \ of the way between centerpoints of 
the links on the unrefined lattice, and should therefore 
be interpolated using 

op') _ qpI v 
oc 2n,n — orj n+N bcl ,n < U n+N bc i-1,& A 

X F 7 ' TI 71 

A J ^n+N bcl -l,h u n+N bcl -l,h ' 

opi) _ Q/Y'ft ttv fro 

°- C/ 2n+l,n — ou 2n,n rj 7i+N bcl ,n u 2n,n 

4- Tf E V U V f (20) 

' u 2n+l,&- Cj n+N b< . 1 + l,n u 2n+l,h > \^ v ) 

where the factor on the LHS is 8 (rather than 4) because 
the lattice E a is a a times the continuum E a , and a v 
is being reduced by a factor of 2. This choice ensures 
that Gauss' law, Eq. (fT7|) . remains identically satisfied at 
even-?) now lattice points (2n,n). 

We define the transverse electric fields emanating from 
the new sites as the average of the parallel transports of 
the electric fields ±1 units of rj away, as shown in Figure 
[3] We parallel transport the indices of E % U l in the same 
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way as we did with the link matrix U l , leading to 
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(21) 



Together with our choice of i?'' fields, this ensures that 
Gauss' law is satisfied also at the new (odd-f/ ncw ) lattice 
sites — up to corrections which arise due to the failure of 
E to commute with the magnetic field, which first arise 
at order a\a^EB 2 . Therefore in the nonabelian context 
Gauss' law will not be identically satisfied at the new in- 
terpolated lattice sites, though the failure is small and 
suppressed by a^a 2 . A final step is needed, in which one 
corrects the electric fields such that Gauss' law becomes 
exact. The optimal choice is to change the electric fields 
by an amount which is strictly the gradient of a scalar 
potential, and with this choice the shift in the electric 
fields, to restore Gauss' law, is unique. We propose to 
use the algorithm for finding this shift, described in de- 
tail in Ref. [ll| . This completes the specification of our 
algorithm for the nonabelian case. 



IV. DISCUSSION 

We have presented a new approach to studying classi- 
cal field theory in a linearly expanding system on the lat- 
tice, with intended applications in the study of early-time 



dynamics after heavy ion collisions. The algorithm allows 
the evolution to proceed to late times in boxes of large 
transverse extent, without encountering the problem that 
the longitudinal (77) lattice spacing becomes coarse. We 
do this by periodically cropping the box in the 77 direction 
and refining the mesh. We have presented a detailed al- 
gorithm both for scalars and for nonabelian gauge fields. 

In practical applications it will probably be necessary 
to start the evolution at very early times, perhaps even 
r ~ a± but in any case r <C L±. We do not think it 
necessary to begin with a range of r\ larger than a few, 
since few excitations will ever propagate over a range of 
77 larger than 1 or 2. Therefore it might make sense to 
begin with a lattice which does have a large hierarchy 
a v <C ax- In this case our cropping and interpolation 
would only begin when a n ~ a±, at times r ~ L±. 

We emphasize again the importance of establishing 
that simulations, particularly of nonabelian gauge fields, 
are really in the large L± limit. This limit is challenging, 
because as we emphasized in the introduction, the screen- 
ing length scale grows as t 1 / 2 . Therefore we expect that 
our approach will actually be necessary to study suffi- 
ciently wide boxes, particularly if the study is to proceed 
to late times. 
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